### Replicate Adams et al with Bootstraps. 
### Produces Results for Table 1 columns 1,2 and 3 (see comments in code on line 123, 136 and 162)
### JBS 06 Sept 2018

set.seed(1234)

rm(list=ls())

### Read in Medians
load('CHES_Bootstrap_1999_2010_medians.Rdata')

final$party_name[final$partyid==103]<-"SP"

parlgov <- read.dta("parlgov_partyids_NEW_shortStata12.dta")
parlgov <- parlgov[,c("partyid","cmp","country_name")]

names(final)
names(parlgov)

nsims <- ncol(final)-5

chesboot.parlgov <- merge(final,parlgov,by.x=c("countryname","partyid"),by.y=c("country_name","partyid"),all.x=TRUE)

chesboot.parlgov<-chesboot.parlgov[complete.cases(chesboot.parlgov),]

chesboot.parlgov$cmp[chesboot.parlgov$cmp==14221] <- 14223
chesboot.parlgov$cmp[chesboot.parlgov$cmp==22220] <- 22210
chesboot.parlgov$cmp[chesboot.parlgov$cmp== 22525|chesboot.parlgov$cmp== 22523] <- 22521
chesboot.parlgov$cmp[chesboot.parlgov$cmp== 51420] <- 51421
chesboot.parlgov$cmp[chesboot.parlgov$cmp== 22524] <- 22320
chesboot.parlgov$cmp[chesboot.parlgov$cmp== 41112] <- 41113
chesboot.parlgov$cmp[chesboot.parlgov$cmp== 31523|chesboot.parlgov$cmp== 31522] <- 31624

chesboot.parlgov <- chesboot.parlgov[chesboot.parlgov$partyid!=308,] # drop CSU

chesboot.parlgov<-chesboot.parlgov[!duplicated(chesboot.parlgov), ]


# Make sure that re-transformed bootstrapped expert scores are used (1-10 scale rather than
# original 1-7 scale, check bootstrap function at top of script)


chesboot.parlgov.for2004<-chesboot.parlgov[chesboot.parlgov$year%in%c(2002,2006),]

interpol2004<-by(chesboot.parlgov.for2004[,c("year",paste("V",1:nsims,sep=""))],list(chesboot.parlgov.for2004$cmp),function(x) {
	if(nrow(x)!=2){
		cat("Check rows: ",rownames(x),"\n")
	}
	else {
			year<-x[,1]
			if (sum(year)==4008){
				dat<-x[,-1]
				colMeans(dat)
			}
		}
	}
)
	
interpol2004<-as.data.frame(do.call("rbind",interpol2004))
interpol2004$year <- 2004
interpol2004$partyid <- rownames(interpol2004)


interpol2009<-by(chesboot.parlgov[chesboot.parlgov$year==2006|chesboot.parlgov$year==2010,c("year",paste("V",1:nsims,sep=""))],list(chesboot.parlgov$cmp[chesboot.parlgov$year==2006|chesboot.parlgov$year==2010]),function (x) {
	if(nrow(x)!=2){
		cat("Check rows: ",rownames(x),"\n")
	}
	else {
	year<-x[,1]
		if(sum(year)==4016){ # check that the two selected years are indeed 2006 and 2010
			dat<-x[,-1]
			return(dat[year==2006,]+3*(dat[year==2010,]-dat[year==2006,])/4)
			}
		}
	}
)

interpol2009<-as.data.frame(do.call("rbind",interpol2009))
interpol2009$year <- 2009
interpol2009$partyid <- rownames(interpol2009)

# merge 1999 CHES data to interpolated 2004 data
# to calculate change

ches1999bootstraps<-chesboot.parlgov[chesboot.parlgov$year==1999,c("cmp",paste("V",1:nsims,sep=""))]
names(ches1999bootstraps)<-c("cmp",paste("1999V",1:nsims,sep=""))

interpol2004.aux<-merge(interpol2004,ches1999bootstraps,by.x="partyid",by.y="cmp")
interpolated_expert.t1.10.boots.2004 <- (interpol2004.aux[,paste("V",1:nsims,sep="")]-interpol2004.aux[,paste("1999V",1:nsims,sep="")])
names(interpolated_expert.t1.10.boots.2004)<-paste("diff",1:nsims,sep="")
rowMeans(interpolated_expert.t1.10.boots.2004)
interpol2004.aux <- cbind(interpol2004.aux,interpolated_expert.t1.10.boots.2004)
interpol2004.aux <- interpol2004.aux[c("partyid",paste("diff",1:nsims,sep=""))]
interpol2004.aux$year<-2004

# merge 2004 interpolated data to interpolated 2009 data
# to calculate change

interpol2004bootstraps<-interpol2004
names(interpol2004bootstraps)<-c(paste("2004V",1:nsims,sep=""),"year","partyid")
interpol2009.aux<-merge(interpol2009,interpol2004bootstraps,by="partyid")
interpolated_expert.t1.10.boots.2009 <- (interpol2009.aux[,paste("V",1:nsims,sep="")]-interpol2009.aux[,paste("2004V",1:nsims,sep="")])
names(interpolated_expert.t1.10.boots.2009)<-paste("diff",1:nsims,sep="")
rowMeans(interpolated_expert.t1.10.boots.2009)
interpol2009.aux <- cbind(interpol2009.aux,interpolated_expert.t1.10.boots.2009)
interpol2009.aux <- interpol2009.aux[c("partyid",paste("diff",1:nsims,sep=""))]
interpol2009.aux$year<-2009

# combine 2004 and 2009 data
interpoldata<-rbind(interpol2004.aux,interpol2009.aux)


# Merge with original replication data file
adams <- read.dta("replication_data_final_AEST.dta")
dim(adams)

adams.boots <- merge(adams,interpoldata,by.x=c("cmpcode","year"),by.y=c("partyid","year"),all.x=T)
dim(adams.boots)

# Replicate Adams et al. 

adams.original.table2.col3 <- lm( party_ch_all_voters_t~ emp_ch1_10+ interpolated_expert_ch1_10,data=adams.boots[is.na(adams.boots$toofew),]) 
# This is the replication of Adams et al without clustering presented in Table 1, column 1 
summary(adams.original.table2.col3)

pm1 <- plm(party_ch_all_voters_t~ emp_ch1_10+ interpolated_expert_ch1_10, data = adams.boots[is.na(adams.boots$toofew),],model= "pooling")

# compute Stata like df-adjustment
G <- length(unique(adams.boots$cmpcode[is.na(adams.boots$toofew)]))
N <- length(adams.boots$cmpcode[is.na(adams.boots$toofew)])
dfa <- (G/(G - 1)) * (N - 1)/pm1$df.residual

# display with cluster VCE and df-adjustment
adams_c_vcov <- dfa * vcovHC(pm1, type = "HC0", cluster = "group", adjust = T)

### These are replicated, clustered se's presented in Table 1, column 2
diag(sqrt(adams_c_vcov))

# Run regressions in each of nsims datasets; report the point predictions in each.
adams.mult.1<-vector('list',nsims)
adams.mult.cf.1<-matrix(rep(NA,nsims*3),nrow=nsims,ncol=3)
adams.mult.se.1<-matrix(rep(NA,nsims*3),nrow=nsims,ncol=3)

# Reduce data to parties used by Adams etal 

adams.boots.red<-adams.boots[is.na(adams.boots$toofew),]
adamsbasedata<-adams.boots.red[,c(1:65)]


for (k in 1:nsims){
	
	currentboot <- cbind(adamsbasedata,interpolated_expert_ch1_10.boot=adams.boots.red[,paste("diff",k,sep="")])	
	
	adams.mult.1[[k]] <- summary(lm(party_ch_all_voters_t~ emp_ch1_10+ interpolated_expert_ch1_10.boot,data=currentboot))
	adams.mult.cf.1[k,] <- as.numeric(adams.mult.1[[k]]$coefficients[,1])
	adams.mult.se.1[k,] <- as.numeric(adams.mult.1[[k]]$coefficients[,2])
	
	currentboot<-NULL
	
	}

### These are the results in Table 1 Column 3 Bootstrapped Median
### Coefficients appear in reverse order of that in the table
### Intercept, Euromanifestos, Experts
mi.meld(adams.mult.cf.1,adams.mult.se.1)




